#include <iostream>
#include <limits>
#include <cmath>
#include <algorithm>
#include <vector>
using namespace std;

class Function
{
public:
    virtual double operator()(double _x) = 0;
};

class Interpolation
{
public:
    virtual double solve(double _x) = 0;
};

class Newton_polynomial : public Function
{
private:
    vector<double> C, X;
    int n;
public:
    Newton_polynomial(vector<double> _C, vector<double> _X) : C(_C), X(_X){}
    double operator()(double _x);
};
double Newton_polynomial::operator()(double _x){
    int i, j;
    n = X.size();
    double p = C[0], pi;
    for (i = 1; i < n; i++){
        pi = C[i];
        for (j = 0; j < i; j++){
            pi = pi * (_x - X[j]);
        }
        p += pi;
    }
    return p;
}

class Newton_formula : public Interpolation
{
private:
    vector<double> X, A;
    Function &f;
    int n;
public:
    Newton_formula(Function &_f, vector<double> _X):f(_f), X(_X){
        n = X.size();
        double t[n][n];
        for (int i = 0; i < n; i++){
            t[i][0] = f(X[i]);
        }
        for (int i = 1; i < n; i++){
            for (int j = i; j < n; j++){
                t[j][i] = (t[j][i - 1] - t[j - 1][i - 1]) / (X[j] - X[j - i]);
            }
        }
        for (int i = 0; i < n; i++){
            A.push_back(t[i][i]);
        }
    }
    double solve(double _x);
    vector<double> Coef();
};

double Newton_formula::solve(double _x){
    Newton_polynomial Pn(A, X);
    return Pn(_x);
}

vector<double> Newton_formula::Coef(){
    return A;
}

class Hermite_polynomial : public Function
{
private:
    vector<double> C, X;
    int n;
public:
    Hermite_polynomial(vector<double> _C, vector<double> _X) : C(_C), X(_X){}
    double operator()(double _x)
    {
        int i, j;
        double p = C[0];
        n = X.size();
        for (i = 1; i < n; i++){
            double pi = C[i];
            for (j = 0; j < i; j++){
                pi = pi * (_x - X[j]);
            }
            p = p + pi;
        }
        return p;
    }
};

class Hermite_formula : public Interpolation
{
private:
    std::vector<double> X, L, A;
    int n;
public:
    Hermite_formula(std::vector<double> _L, std::vector<double> _X): L(_L), X(_X){
        n = X.size();
        int m = 1;
        double t[n][n] = {{0}};
        for (int i = 0; i < n; i++){
            for (int j = 0; j < m; j++){
                t[i][j] = X[i + j + 1 - m];
            }
            if (L[i] == L[i + 1]){
                m = m + 1;
            }
            else{
                m = 1;
            }
        }
        for (int i = 1; i < n; i++){
            for (int j = i; j < n; j++){
                if (t[j][i] == 0)
                    t[j][i] = (t[j][i - 1] - t[j - 1][i - 1]) / (L[j] - L[j - i]);
            }
        }
        for (int i = 0; i < n; i++){
            A.push_back(t[i][i]);
        }
    }
    double solve(double _x);
    double diff(double _x);
    vector<double> Coef();
};

double Hermite_formula::solve(double _x){
    Newton_polynomial Pn(A, L);
    return Pn(_x);
}

double Hermite_formula::diff(double _x){
    double epsilon = 0.000001;
    Newton_polynomial Pn(A, L);
    return (Pn(_x + epsilon) - Pn(_x - epsilon)) / (2 * epsilon);
}

vector<double> Hermite_formula::Coef(){
    return A;
}
